#ifndef PCPS_BLAS_LAPACK_HEADER
#define PCPS_BLAS_LAPACK_HEADER

#include <complex>
#include <typeinfo>

#include <boost/format.hpp>
#include <src/exception.h>

// *** Custom BLAS-like functions ***

namespace pcps {

  // complex conjugate of a vector
  inline void xconj(const int n, double * const y, const int incy) {
  }
  inline void xconj(const int n, std::complex<double> * const y, const int incy) {
    const int n5 = n - n % 5;
    int i = 0;
    if (incy == 1) {
      for ( ; i < n5; i += 5) {
        y[i+0] = std::conj(y[i+0]);
        y[i+1] = std::conj(y[i+1]);
        y[i+2] = std::conj(y[i+2]);
        y[i+3] = std::conj(y[i+3]);
        y[i+4] = std::conj(y[i+4]);
      }
      for ( ; i < n; i += 1)
        y[i+0] = std::conj(y[i+0]);
    } else {
      for ( ; i < n5; i += 5) {
        y[(i+0)*incy] = std::conj(y[(i+0)*incy]);
        y[(i+1)*incy] = std::conj(y[(i+1)*incy]);
        y[(i+2)*incy] = std::conj(y[(i+2)*incy]);
        y[(i+3)*incy] = std::conj(y[(i+3)*incy]);
        y[(i+4)*incy] = std::conj(y[(i+4)*incy]);
      }
      for ( ; i < n; i += 1)
        y[(i+0)*incy] = std::conj(y[(i+0)*incy]);
    }
  }

}

// *** MKL Interface ***

#ifdef HAVE_MKL

#include <mkl.h>

namespace pcps {

  // multiply vector by a constant
  inline void dscal(const MKL_INT n, const double a, double *x, const MKL_INT incx) {
    DSCAL(&n, &a, x, &incx);
  }
  inline void xscal(const MKL_INT n, const double a, double *x, const MKL_INT incx) {
    DSCAL(&n, &a, x, &incx);
  }
  inline void zscal(const MKL_INT n, const std::complex<double> a, std::complex<double> *x, const MKL_INT incx) {
    ZSCAL(&n, (const MKL_Complex16 *)&a, (MKL_Complex16 *)x, &incx);
  }
  inline void xscal(const MKL_INT n, const std::complex<double> a, std::complex<double> *x, const MKL_INT incx) {
    ZSCAL(&n, (const MKL_Complex16 *)&a, (MKL_Complex16 *)x, &incx);
  }

  // copy a vector
  inline void dcopy(const MKL_INT n, const double *x, const MKL_INT incx, double *y, MKL_INT incy) {
    DCOPY(&n, x, &incx, y, &incy);
  }
  inline void xcopy(const MKL_INT n, const double *x, const MKL_INT incx, double *y, MKL_INT incy) {
    DCOPY(&n, x, &incx, y, &incy);
  }
  inline void zcopy(const MKL_INT n, const std::complex<double> *x, const MKL_INT incx, std::complex<double> *y, MKL_INT incy) {
    ZCOPY(&n, (const MKL_Complex16 *)x, &incx, (MKL_Complex16 *)y, &incy);
  }
  inline void xcopy(const MKL_INT n, const std::complex<double> *x, const MKL_INT incx, std::complex<double> *y, MKL_INT incy) {
    ZCOPY(&n, (const MKL_Complex16 *)x, &incx, (MKL_Complex16 *)y, &incy);
  }

  // vector addition
  inline void daxpy(const MKL_INT n, const double a, const double *x, const MKL_INT incx, double *y, const MKL_INT incy) {
    DAXPY(&n, &a, x, &incx, y, &incy);
  }
  inline void xaxpy(const MKL_INT n, const double a, const double *x, const MKL_INT incx, double *y, const MKL_INT incy) {
    DAXPY(&n, &a, x, &incx, y, &incy);
  }
  inline void zaxpy(const MKL_INT n, const std::complex<double> a, const std::complex<double> *x, const MKL_INT incx,
                    std::complex<double> *y, const MKL_INT incy) {
    ZAXPY(&n, (const MKL_Complex16 *)&a, (const MKL_Complex16 *)x, &incx, (MKL_Complex16 *)y, &incy);
  }
  inline void xaxpy(const MKL_INT n, const std::complex<double> a, const std::complex<double> *x, const MKL_INT incx,
                    std::complex<double> *y, const MKL_INT incy) {
    ZAXPY(&n, (const MKL_Complex16 *)&a, (const MKL_Complex16 *)x, &incx, (MKL_Complex16 *)y, &incy);
  }

  // dot product
  inline double ddot(const MKL_INT n, const double *x, const MKL_INT incx, const double *y, const MKL_INT incy) {
    return DDOT(&n, x, &incx, y, &incy);
  }
  inline double xdot(const MKL_INT n, const double *x, const MKL_INT incx, const double *y, const MKL_INT incy) {
    return DDOT(&n, x, &incx, y, &incy);
  }
  inline std::complex<double> zdot(const MKL_INT n, const std::complex<double> *x, const MKL_INT incx, const std::complex<double> *y, const MKL_INT incy) {
    std::complex<double> retval;
    ZDOTU((MKL_Complex16 *)&retval, &n, (const MKL_Complex16 *)x, &incx, (const MKL_Complex16 *)y, &incy);
    return retval;
  }
  inline std::complex<double> xdot(const MKL_INT n, const std::complex<double> *x, const MKL_INT incx, const std::complex<double> *y, const MKL_INT incy) {
    std::complex<double> retval;
    ZDOTU((MKL_Complex16 *)&retval, &n, (const MKL_Complex16 *)x, &incx, (const MKL_Complex16 *)y, &incy);
    return retval;
  }

  // dot product with first vector conjugated
  inline double ddotc(const MKL_INT n, const double *x, const MKL_INT incx, const double *y, const MKL_INT incy) {
    return DDOT(&n, x, &incx, y, &incy);
  }
  inline double xdotc(const MKL_INT n, const double *x, const MKL_INT incx, const double *y, const MKL_INT incy) {
    return DDOT(&n, x, &incx, y, &incy);
  }
  inline std::complex<double> zdotc(const MKL_INT n, const std::complex<double> *x, const MKL_INT incx, const std::complex<double> *y, const MKL_INT incy) {
    std::complex<double> retval;
    ZDOTC((MKL_Complex16 *)&retval, &n, (const MKL_Complex16 *)x, &incx, (const MKL_Complex16 *)y, &incy);
    return retval;
  }
  inline std::complex<double> xdotc(const MKL_INT n, const std::complex<double> *x, const MKL_INT incx, const std::complex<double> *y, const MKL_INT incy) {
    std::complex<double> retval;
    ZDOTC((MKL_Complex16 *)&retval, &n, (const MKL_Complex16 *)x, &incx, (const MKL_Complex16 *)y, &incy);
    return retval;
  }

  // matrix multiplication
  inline void dgemm(const char transa, const char transb, const MKL_INT m, const MKL_INT n, const MKL_INT k,
                    const double alpha, const double *a, const MKL_INT lda, const double *b, const MKL_INT ldb,
                    const double beta, double *c, const MKL_INT ldc) {
    DGEMM(&transa, &transb, &m, &n, &k, &alpha, a, &lda, b, &ldb, &beta, c, &ldc);
  }
  inline void xgemm(const char transa, const char transb, const MKL_INT m, const MKL_INT n, const MKL_INT k,
                    const double alpha, const double *a, const MKL_INT lda, const double *b, const MKL_INT ldb,
                    const double beta, double *c, const MKL_INT ldc) {
    DGEMM(&transa, &transb, &m, &n, &k, &alpha, a, &lda, b, &ldb, &beta, c, &ldc);
  }
  inline void zgemm(const char transa, const char transb, const MKL_INT m, const MKL_INT n, const MKL_INT k,
                    const std::complex<double> alpha, const std::complex<double> *a, const MKL_INT lda,
                    const std::complex<double> *b, const MKL_INT ldb,
                    const std::complex<double> beta, std::complex<double> *c, const MKL_INT ldc) {
    ZGEMM(&transa, &transb, &m, &n, &k, (const MKL_Complex16 *)&alpha, (const MKL_Complex16 *)a, &lda,
          (const MKL_Complex16 *)b, &ldb, (const MKL_Complex16 *)&beta, (MKL_Complex16 *)c, &ldc);
  }
  inline void xgemm(const char transa, const char transb, const MKL_INT m, const MKL_INT n, const MKL_INT k,
                    const std::complex<double> alpha, const std::complex<double> *a, const MKL_INT lda,
                    const std::complex<double> *b, const MKL_INT ldb,
                    const std::complex<double> beta, std::complex<double> *c, const MKL_INT ldc) {
    ZGEMM(&transa, &transb, &m, &n, &k, (const MKL_Complex16 *)&alpha, (const MKL_Complex16 *)a, &lda,
          (const MKL_Complex16 *)b, &ldb, (const MKL_Complex16 *)&beta, (MKL_Complex16 *)c, &ldc);
  }

  // solve linear equations
  inline void dgesv(MKL_INT n, MKL_INT nrhs, double *a, MKL_INT lda,
                    MKL_INT *ipiv, double *b, MKL_INT ldb, MKL_INT & info) {
    DGESV(&n, &nrhs, a, &lda, ipiv, b, &ldb, &info);
  }
  inline void xgesv(MKL_INT n, MKL_INT nrhs, double *a, MKL_INT lda,
                    MKL_INT *ipiv, double *b, MKL_INT ldb, MKL_INT & info) {
    pcps::dgesv(n, nrhs, a, lda, ipiv, b, ldb, info);
  }
  inline void zgesv(MKL_INT n, MKL_INT nrhs, std::complex<double> *a, MKL_INT lda,
                    MKL_INT *ipiv, std::complex<double> *b, MKL_INT ldb, MKL_INT & info) {
    ZGESV(&n, &nrhs, (MKL_Complex16 *)a, &lda, ipiv, (MKL_Complex16 *)b, &ldb, &info);
  }
  inline void xgesv(MKL_INT n, MKL_INT nrhs, std::complex<double> *a, MKL_INT lda,
                    MKL_INT *ipiv, std::complex<double> *b, MKL_INT ldb, MKL_INT & info) {
    pcps::zgesv(n, nrhs, a, lda, ipiv, b, ldb, info);
  }

  // solve symmetric/hermitian eigensystem
  inline void dsyev(char jobz, char uplo, MKL_INT n, double *a, MKL_INT lda,
                    double * w, double * work, MKL_INT lwork, MKL_INT & info) {
    DSYEV(&jobz, &uplo, &n, a, &lda, w, work, &lwork, &info);
  }
  inline void xsyev(char jobz, char uplo, MKL_INT n, double *a, MKL_INT lda,
                    double * w, double * work, MKL_INT lwork, double * rwork, MKL_INT & info) {
    pcps::dsyev(jobz, uplo, n, a, lda, w, work, lwork, info);
  }
  inline void zheev(char jobz, char uplo, MKL_INT n, std::complex<double> *a, MKL_INT lda,
                    double * w, std::complex<double> * work, MKL_INT lwork, double * rwork, MKL_INT & info) {
    ZHEEV(&jobz, &uplo, &n, (MKL_Complex16 *)a, &lda, w, (MKL_Complex16 *)work, &lwork, rwork, &info);
  }
  inline void xsyev(char jobz, char uplo, MKL_INT n, std::complex<double> *a, MKL_INT lda,
                    double * w, std::complex<double> * work, MKL_INT lwork, double * rwork, MKL_INT & info) {
    pcps::zheev(jobz, uplo, n, a, lda, w, work, lwork, rwork, info);
  }

  // solve non-symmetric/non-hermitian eigensystem
  inline void dgeev(char jobvl, char jobvr, MKL_INT n, double *a, MKL_INT lda,
                    double *wr, double *wi, double *vl, MKL_INT ldvl, double *vr, MKL_INT ldvr,
                    double *work, MKL_INT lwork, MKL_INT & info) {
    DGEEV(&jobvl, &jobvr, &n, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr, work, &lwork, &info);
  }
  inline void xgeev(char jobvl, char jobvr, MKL_INT n, double *a, MKL_INT lda,
                    double *w, double *wr, double *wi, double *vl, MKL_INT ldvl, double *vr, MKL_INT ldvr,
                    double *work, MKL_INT lwork, double *rwork, MKL_INT & info) {
    pcps::dgeev(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info);
  }
  inline void zgeev(char jobvl, char jobvr, MKL_INT n, std::complex<double> *a, MKL_INT lda,
                    std::complex<double> *w, std::complex<double> *vl, MKL_INT ldvl, std::complex<double> *vr, MKL_INT ldvr,
                    std::complex<double> *work, MKL_INT lwork, double *rwork, MKL_INT & info) {
    ZGEEV(&jobvl, &jobvr, &n, (MKL_Complex16 *)a, &lda, (MKL_Complex16 *)w,
          (MKL_Complex16 *)vl, &ldvl, (MKL_Complex16 *)vr, &ldvr, (MKL_Complex16 *)work, &lwork, rwork, &info);
  }
  inline void xgeev(char jobvl, char jobvr, MKL_INT n, std::complex<double> *a, MKL_INT lda,
                    std::complex<double> *w, double *wr, double *wi, std::complex<double> *vl, MKL_INT ldvl, std::complex<double> *vr, MKL_INT ldvr,
                    std::complex<double> *work, MKL_INT lwork, double *rwork, MKL_INT & info) {
    pcps::zgeev(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info);
    for (int i = 0; i < n; i++) {
      wr[i] = w[i].real();
      wi[i] = w[i].imag();
    }
  }

  // singular value decomposition
  inline void dgesvd(char jobu, char jobvt, MKL_INT m, MKL_INT n, double *a, MKL_INT lda,
                     double *s, double *u, MKL_INT ldu, double *vt, MKL_INT ldvt,
                     double *work, MKL_INT lwork, MKL_INT & info) {
    DGESVD(&jobu, &jobvt, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, &info);
  }
  inline void xgesvd(char jobu, char jobvt, MKL_INT m, MKL_INT n, double *a, MKL_INT lda,
                     double *s, double *u, MKL_INT ldu, double *vt, MKL_INT ldvt,
                     double *work, MKL_INT lwork, double *rwork, MKL_INT & info) {
    pcps::dgesvd(jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, info);
  }
  inline void zgesvd(char jobu, char jobvt, MKL_INT m, MKL_INT n, std::complex<double> *a, MKL_INT lda,
                     double *s, std::complex<double> *u, MKL_INT ldu, std::complex<double> *vt, MKL_INT ldvt,
                     std::complex<double> *work, MKL_INT lwork, double *rwork, MKL_INT & info) {
    ZGESVD(&jobu, &jobvt, &m, &n, (MKL_Complex16 *)a, &lda, s, (MKL_Complex16 *)u, &ldu, (MKL_Complex16 *)vt, &ldvt,
           (MKL_Complex16 *)work, &lwork, rwork, &info);
  }
  inline void xgesvd(char jobu, char jobvt, MKL_INT m, MKL_INT n, std::complex<double> *a, MKL_INT lda,
                     double *s, std::complex<double> *u, MKL_INT ldu, std::complex<double> *vt, MKL_INT ldvt,
                     std::complex<double> *work, MKL_INT lwork, double *rwork, MKL_INT & info) {
    pcps::zgesvd(jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, rwork, info);
  }

  // LU decomposition
  inline void dgetrf(MKL_INT m, MKL_INT n, double *a, MKL_INT lda, MKL_INT * ipiv, MKL_INT & info) {
    DGETRF(&m, &n, a, &lda, ipiv, &info);
  }
  inline void xgetrf(MKL_INT m, MKL_INT n, double *a, MKL_INT lda, MKL_INT * ipiv, MKL_INT & info) {
    pcps::dgetrf(m, n, a, lda, ipiv, info);
  }
  inline void zgetrf(MKL_INT m, MKL_INT n, std::complex<double> *a, MKL_INT lda, MKL_INT * ipiv, MKL_INT & info) {
    ZGETRF(&m, &n, (MKL_Complex16 *)a, &lda, ipiv, &info);
  }
  inline void xgetrf(MKL_INT m, MKL_INT n, std::complex<double> *a, MKL_INT lda, MKL_INT * ipiv, MKL_INT & info) {
    pcps::zgetrf(m, n, a, lda, ipiv, info);
  }

}

#endif

// *** Standard BLAS/LAPACK Interface ***

#ifndef HAVE_MKL

#define FFUNC(a) a ## _

// declare functions we need from the BLAS/LAPACK library
extern "C" {
  void FFUNC(dscal)(const int *n, const double *a, double *x, const int *incx);
  void FFUNC(zscal)(const int *n, const std::complex<double> *a, std::complex<double> *x, const int *incx);
  void FFUNC(dcopy)(const int *n, const double * x, const int *incx, double *y, const int *incy);
  void FFUNC(zcopy)(const int *n, const std::complex<double> * x, const int *incx, std::complex<double> *y, const int *incy);
  void FFUNC(daxpy)(const int *n, const double *a, const double *x, const int *incx, double *y, const int *incy);
  void FFUNC(zaxpy)(const int *n, const std::complex<double> *a, const std::complex<double> *x, const int *incx,
                                                                       std::complex<double> *y, const int *incy);
  double FFUNC(ddot)(const int *n, const double *x, const int *incx, const double *y, const int *incy);
  std::complex<double> FFUNC(zdotu)(const int *n, const std::complex<double> *x, const int *incx,
                                                  const std::complex<double> *y, const int *incy);
  std::complex<double> FFUNC(zdotc)(const int *n, const std::complex<double> *x, const int *incx,
                                                  const std::complex<double> *y, const int *incy);
  void FFUNC(dgemm)(const char *transa, const char *transb, const int *m, const int *n, const int *k,
                    const double *alpha, const double *a, const int *lda, const double *b, const int *ldb,
                    const double *beta, double *c, const int *ldc);
  void FFUNC(zgemm)(const char *transa, const char *transb, const int *m, const int *n, const int *k,
                    const std::complex<double> *alpha,
                    const std::complex<double> *a, const int *lda,
                    const std::complex<double> *b, const int *ldb,
                    const std::complex<double> *beta,
                          std::complex<double> *c, const int *ldc);
  void FFUNC(dgesv)(const int *n, const int *nrhs, double *a, const int *lda, int *ipiv,
                    double *b, const int *ldb, int *info);
  void FFUNC(zgesv)(const int *n, const int *nrhs, std::complex<double> *a, const int *lda, int *ipiv,
                    std::complex<double> *b, const int *ldb, int *info);
  void FFUNC(dsyev)(char *jobz, char *uplo, int *n, double *a, int *lda, double *w,
                    double *work, int *lwork, int *info);
  void FFUNC(zheev)(char *jobz, char *uplo, int *n, std::complex<double> *a, int *lda, double *w,
                    std::complex<double> *work, int *lwork, double *rwork, int *info);
  void FFUNC(dgeev)(char *jobvl, char *jobvr, int *n, double *a, int *lda, double *wr, double *wi,
                    double *vl, int *ldvl, double *vr, int *ldvr, double *work, int *lwork, int *info);
  void FFUNC(zgeev)(char *jobvl, char *jobvr, int *n, std::complex<double> *a, int *lda, std::complex<double> *w,
                    std::complex<double> *vl, int *ldvl, std::complex<double> *vr, int *ldvr,
                    std::complex<double> *work, int *lwork, double *rwork, int *info);
  void FFUNC(dgesvd)(char *jobu, char *jobvt, int *m, int *n, double *a, int *lda, double *s,
                     double *u, int *ldu, double *vt, int *ldvt, double *work, int *lwork, int *info);
  void FFUNC(zgesvd)(char *jobu, char *jobvt, int *m, int *n, std::complex<double> *a, int *lda, double *s,
                     std::complex<double> *u, int *ldu, std::complex<double> *vt, int *ldvt,
                     std::complex<double> *work, int *lwork, double *rwork, int *info);
  void FFUNC(dgetrf)(int *m, int *n, double *a, int *lda, int *ipiv, int *info);
  void FFUNC(zgetrf)(int *m, int *n, std::complex<double> *a, int *lda, int *ipiv, int *info);
}

// define interface functions
namespace pcps {

  // multiply vector by a constant
  inline void dscal(const int n, const double a, double *x, const int incx) {
    FFUNC(dscal)(&n, &a, x, &incx);
  }
  inline void xscal(const int n, const double a, double *x, const int incx) {
    FFUNC(dscal)(&n, &a, x, &incx);
  }
  inline void zscal(const int n, const std::complex<double> a, std::complex<double> *x, const int incx) {
    FFUNC(zscal)(&n, &a, x, &incx);
  }
  inline void xscal(const int n, const std::complex<double> a, std::complex<double> *x, const int incx) {
    FFUNC(zscal)(&n, &a, x, &incx);
  }

  // copy a vector
  inline void dcopy(const int n, const double *x, const int incx, double *y, int incy) {
    FFUNC(dcopy)(&n, x, &incx, y, &incy);
  }
  inline void xcopy(const int n, const double *x, const int incx, double *y, int incy) {
    FFUNC(dcopy)(&n, x, &incx, y, &incy);
  }
  inline void zcopy(const int n, const std::complex<double> *x, const int incx, std::complex<double> *y, int incy) {
    FFUNC(zcopy)(&n, x, &incx, y, &incy);
  }
  inline void xcopy(const int n, const std::complex<double> *x, const int incx, std::complex<double> *y, int incy) {
    FFUNC(zcopy)(&n, x, &incx, y, &incy);
  }

  // vector addition
  inline void daxpy(const int n, const double a, const double *x, const int incx, double *y, const int incy) {
    FFUNC(daxpy)(&n, &a, x, &incx, y, &incy);
  }
  inline void xaxpy(const int n, const double a, const double *x, const int incx, double *y, const int incy) {
    FFUNC(daxpy)(&n, &a, x, &incx, y, &incy);
  }
  inline void zaxpy(const int n, const std::complex<double> a, const std::complex<double> *x, const int incx,
                    std::complex<double> *y, const int incy) {
    FFUNC(zaxpy)(&n, &a, x, &incx, y, &incy);
  }
  inline void xaxpy(const int n, const std::complex<double> a, const std::complex<double> *x, const int incx,
                    std::complex<double> *y, const int incy) {
    FFUNC(zaxpy)(&n, &a, x, &incx, y, &incy);
  }

  // dot product
  inline double ddot(const int n, const double *x, const int incx, const double *y, const int incy) {
    return FFUNC(ddot)(&n, x, &incx, y, &incy);
  }
  inline double xdot(const int n, const double *x, const int incx, const double *y, const int incy) {
    return FFUNC(ddot)(&n, x, &incx, y, &incy);
  }
  inline std::complex<double> zdot(const int n, const std::complex<double> *x, const int incx, const std::complex<double> *y, const int incy) {
    return FFUNC(zdotu)(&n, x, &incx, y, &incy);
  }
  inline std::complex<double> xdot(const int n, const std::complex<double> *x, const int incx, const std::complex<double> *y, const int incy) {
    return FFUNC(zdotu)(&n, x, &incx, y, &incy);
  }

  // dot product with first vector conjugated
  inline double ddotc(const int n, const double *x, const int incx, const double *y, const int incy) {
    return FFUNC(ddot)(&n, x, &incx, y, &incy);
  }
  inline double xdotc(const int n, const double *x, const int incx, const double *y, const int incy) {
    return FFUNC(ddot)(&n, x, &incx, y, &incy);
  }
  inline std::complex<double> zdotc(const int n, const std::complex<double> *x, const int incx, const std::complex<double> *y, const int incy) {
    return FFUNC(zdotc)(&n, x, &incx, y, &incy);
  }
  inline std::complex<double> xdotc(const int n, const std::complex<double> *x, const int incx, const std::complex<double> *y, const int incy) {
    return FFUNC(zdotc)(&n, x, &incx, y, &incy);
  }

  // matrix multiplication
  inline void dgemm(const char transa, const char transb, const int m, const int n, const int k,
                    const double alpha, const double *a, const int lda, const double *b, const int ldb,
                    const double beta, double *c, const int ldc) {
    FFUNC(dgemm)(&transa, &transb, &m, &n, &k, &alpha, a, &lda, b, &ldb, &beta, c, &ldc);
  }
  inline void xgemm(const char transa, const char transb, const int m, const int n, const int k,
                    const double alpha, const double *a, const int lda, const double *b, const int ldb,
                    const double beta, double *c, const int ldc) {
    FFUNC(dgemm)(&transa, &transb, &m, &n, &k, &alpha, a, &lda, b, &ldb, &beta, c, &ldc);
  }
  inline void zgemm(const char transa, const char transb, const int m, const int n, const int k,
                    const std::complex<double> alpha, const std::complex<double> *a, const int lda,
                    const std::complex<double> *b, const int ldb,
                    const std::complex<double> beta, std::complex<double> *c, const int ldc) {
    FFUNC(zgemm)(&transa, &transb, &m, &n, &k, &alpha, a, &lda, b, &ldb, &beta, c, &ldc);
  }
  inline void xgemm(const char transa, const char transb, const int m, const int n, const int k,
                    const std::complex<double> alpha, const std::complex<double> *a, const int lda,
                    const std::complex<double> *b, const int ldb,
                    const std::complex<double> beta, std::complex<double> *c, const int ldc) {
    FFUNC(zgemm)(&transa, &transb, &m, &n, &k, &alpha, a, &lda, b, &ldb, &beta, c, &ldc);
  }

  // solve linear equations
  inline void dgesv(int n, int nrhs, double *a, int lda,
                    int *ipiv, double *b, int ldb, int & info) {
    FFUNC(dgesv)(&n, &nrhs, a, &lda, ipiv, b, &ldb, &info);
  }
  inline void xgesv(int n, int nrhs, double *a, int lda,
                    int *ipiv, double *b, int ldb, int & info) {
    pcps::dgesv(n, nrhs, a, lda, ipiv, b, ldb, info);
  }
  inline void zgesv(int n, int nrhs, std::complex<double> *a, int lda,
                    int *ipiv, std::complex<double> *b, int ldb, int & info) {
    FFUNC(zgesv)(&n, &nrhs, a, &lda, ipiv, b, &ldb, &info);
  }
  inline void xgesv(int n, int nrhs, std::complex<double> *a, int lda,
                    int *ipiv, std::complex<double> *b, int ldb, int & info) {
    pcps::zgesv(n, nrhs, a, lda, ipiv, b, ldb, info);
  }

  // solve symmetric/hermitian eigensystem
  inline void dsyev(char jobz, char uplo, int n, double *a, int lda,
                    double * w, double * work, int lwork, int & info) {
    FFUNC(dsyev)(&jobz, &uplo, &n, a, &lda, w, work, &lwork, &info);
  }
  inline void xsyev(char jobz, char uplo, int n, double *a, int lda,
                    double * w, double * work, int lwork, double * rwork, int & info) {
    pcps::dsyev(jobz, uplo, n, a, lda, w, work, lwork, info);
  }
  inline void zheev(char jobz, char uplo, int n, std::complex<double> *a, int lda,
                    double * w, std::complex<double> * work, int lwork, double * rwork, int & info) {
    FFUNC(zheev)(&jobz, &uplo, &n, a, &lda, w, work, &lwork, rwork, &info);
  }
  inline void xsyev(char jobz, char uplo, int n, std::complex<double> *a, int lda,
                    double * w, std::complex<double> * work, int lwork, double * rwork, int & info) {
    pcps::zheev(jobz, uplo, n, a, lda, w, work, lwork, rwork, info);
  }

  // solve non-symmetric/non-hermitian eigensystem
  inline void dgeev(char jobvl, char jobvr, int n, double *a, int lda,
                    double *wr, double *wi, double *vl, int ldvl, double *vr, int ldvr,
                    double *work, int lwork, int & info) {
    FFUNC(dgeev)(&jobvl, &jobvr, &n, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr, work, &lwork, &info);
  }
  inline void xgeev(char jobvl, char jobvr, int n, double *a, int lda,
                    double *w, double *wr, double *wi, double *vl, int ldvl, double *vr, int ldvr,
                    double *work, int lwork, double *rwork, int & info) {
    pcps::dgeev(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info);
  }
  inline void zgeev(char jobvl, char jobvr, int n, std::complex<double> *a, int lda,
                    std::complex<double> *w, std::complex<double> *vl, int ldvl, std::complex<double> *vr, int ldvr,
                    std::complex<double> *work, int lwork, double *rwork, int & info) {
    FFUNC(zgeev)(&jobvl, &jobvr, &n, a, &lda, w, vl, &ldvl, vr, &ldvr, work, &lwork, rwork, &info);
  }
  inline void xgeev(char jobvl, char jobvr, int n, std::complex<double> *a, int lda,
                    std::complex<double> *w, double *wr, double *wi, std::complex<double> *vl, int ldvl,
                    std::complex<double> *vr, int ldvr, std::complex<double> *work, int lwork,
                    double *rwork, int & info) {
    pcps::zgeev(jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr, work, lwork, rwork, info);
    for (int i = 0; i < n; i++) {
      wr[i] = w[i].real();
      wi[i] = w[i].imag();
    }
  }

  // singular value decomposition
  inline void dgesvd(char jobu, char jobvt, int m, int n, double *a, int lda,
                     double *s, double *u, int ldu, double *vt, int ldvt,
                     double *work, int lwork, int & info) {
    FFUNC(dgesvd)(&jobu, &jobvt, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, &info);
  }
  inline void xgesvd(char jobu, char jobvt, int m, int n, double *a, int lda,
                     double *s, double *u, int ldu, double *vt, int ldvt,
                     double *work, int lwork, double *rwork, int & info) {
    pcps::dgesvd(jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, info);
  }
  inline void zgesvd(char jobu, char jobvt, int m, int n, std::complex<double> *a, int lda,
                     double *s, std::complex<double> *u, int ldu, std::complex<double> *vt, int ldvt,
                     std::complex<double> *work, int lwork, double *rwork, int & info) {
    FFUNC(zgesvd)(&jobu, &jobvt, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, rwork, &info);
  }
  inline void xgesvd(char jobu, char jobvt, int m, int n, std::complex<double> *a, int lda,
                     double *s, std::complex<double> *u, int ldu, std::complex<double> *vt, int ldvt,
                     std::complex<double> *work, int lwork, double *rwork, int & info) {
    pcps::zgesvd(jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, rwork, info);
  }

  // LU decomposition
  inline void dgetrf(int m, int n, double *a, int lda, int * ipiv, int & info) {
    FFUNC(dgetrf)(&m, &n, a, &lda, ipiv, &info);
  }
  inline void xgetrf(int m, int n, double *a, int lda, int * ipiv, int & info) {
    pcps::dgetrf(m, n, a, lda, ipiv, info);
  }
  inline void zgetrf(int m, int n, std::complex<double> *a, int lda, int * ipiv, int & info) {
    FFUNC(zgetrf)(&m, &n, a, &lda, ipiv, &info);
  }
  inline void xgetrf(int m, int n, std::complex<double> *a, int lda, int * ipiv, int & info) {
    pcps::zgetrf(m, n, a, lda, ipiv, info);
  }

}

#undef FFUNC

#endif

#endif
